PCT/JP03/13132 

26.1 2.03 



PA 1097492 



WiPO 



pen 



UNITED STATES DEPARTMENT OF COMMERCE 
United States Patent and Trademark Office 

November 24, 2003 

TfflS IS TO CERTIFY THAT ANNEXED HERETO IS A TRUE COPY FROM 
THE RECORDS OF THE UNITED STATES PATENT AND TRADEMARK 
OFFICE OF THOSE PAPERS OF THE BELOW IDENTIFIED PATENT 
APPLICATION THAT MET THE REQUIREMENTS TO BE GRANTED A 
FILING DATE UNDER 35 USC 111. 

APPLICATION NUMBER: 60/418,313 
FILING DATE: October 15, 2002 




PRIORITY DOCUMENT 

SUBMITTED OR TRANSMITTED tm 
COMPLIANCE WITH 
RULE 17.1(a) OR (b) 



By Autliority of the 
OMMISSIONER OF PATENTS AND TRADEMARKS 

p. SVJ55UN 
Certifying Officer 



c=>l 



o 

: ^ 
!c 





Pleasa typo a phis sign {*) tnside this box • 



PTO/SB/16 (02-01) 
Approved for use throughl 0/31/2002 OMB 0651-0032 
U S. Patent and Trademark Office. U S DEPARTMEMT OF COMM£RCI|^^i ^^ 



Under the Papemork Reduction Act of 1995. no persons are requred to respond to a icOecSion of infonnation unless il displays a valid OMB control num^fJJ^g 



PROVISIONAL APPLICATION FOR PATENT COVER SHEET 

This Is a request for filing a PROVISIONAL APPUCATION FOR PATEN T under 37 CFR 1.53(c). 
I Express Maill^brfN^ — ^—i ^^^^^^^^dL 




Ghren Name (first and mlddlB Qf anW) 



Akinori 



iriVENTOR(S) 



Family Name or Sumanrte 



Nishihara 



ReskJencB 
(City and either State or Foreign (kiuntiy) 



374.4.505,lmainakamachi.Nakahara-ku, 
Kawasaki-shLKanagawaJapan 



13 Aaditsonalaiventora are being named oniha separately numbered sheets attached hereto 

TITLE OF THE INVENTION (280 charac ters maxi _ 

Algorithm for Exacl Computation of Coefficients of Universal MaximaOy Flat FIR Filters 



Ofrecf aO conespondenca fo: 
\ I Customer Number 
OR 



CORRESPONDENCE ADDRESS 
► 



Type Customer Number here 



Piece Customer Number 
Bar Code Label hare 



^ FUmOr 



Individual Name 



Address 



Akinbri Nishihara 



374>4-5bSJmainakarrtachi,Nakahara:ku,KaA^^^ 



1 Addfess 
1 CHy 


Kawasaki-shi 


state 


Kanagawa 


ZIP 


211-0065 1 




Jaoan 


Teleohone 


+81^3-5734-3232 




+81-3-5734-2994 ' 1 



Specification Number of Pages | T6 1 
1^ Drawing(s) Number (^Sheets Q 
I I Application Data Sheet See 37 CFR 1 76 



ENCLOSED APPUCATION PARTS (c heck all that agp/y) 



3 



\ \ (XKs). Number 
I \ Other (spedW 



METHOD 



OF PAYMENT OF FILING FEES FOR THIS PROVISIONAL APPLICATION FOR PATEhfr 



Applicant claims smaU entity status See 37 CFR 1 ^7 
A check or money order is enclosed to cover the fiflng fees 
The Comnrilssioner is hereby authonzed to charge filing r 
fees or credit any overpayment to Deposit Account Number L 
Payment by credit card Form PTO-2038 is attached. 



FILING FEE 
AMOUNT (S^ 



$80 



The invention was mads by an agency of the United States Government or under a contract with an agency of the 
United Stales GovemmenL 
IE^ No 

D Yes. the name of the us Governmenl egancy end the Government contract number are _ 



Respectfully 

signaturI 




Date 



\0i 4 I 'CC. 



TVPED or PRINTED NAME . 



Akinori Nishihara 



TELEPHONE 



+81-3-5734-3232 



REGISTRATION NO. 
Qf appropriate) 
Docket Number 



USE ONLY FOR FILING A PROVISIONAL APPLICATION FOR PATENT 

This conecaon of WbimaOon Is required by 37 CFR 1.51 The WoimaUon b iged by thepubllc to ^^^^S?J^?Jil ff^l5^\n 
oSlfftf^Xl S2S ^SSfmSfTtS^n fte amount ^ ti^e you reffilre lb complete andto suggs^^ 

Comm'lssloner for Patents. Washlniston, D.C. 2Q231. 



CoDV provided bv USPTO from the PACR Imaae Database on 11/21/2003 



PROVISIONAL APPLICATION COVER SHEET 

Additional Page 



Under Iha Papenrork Reduction Act a> 188S, no poisons ara reqi 



PToraayis (02-01) 

Approved for use through 10im«002 OMB 0651-0032 
U S PaterU and Trademark Office, U S DEPARTMEWT OF COMMERCE 
torespondtoa^Ue^ion^iT^or^ 



Dodcet Number 



Type a plus sign (■♦') 
inside Ihis box 



Given Name (first and middle Df anyP 



INVENTOR(Sl/APPUCANT(8) 



Family or Surname 



Residence 
f City and either State or Foreign Countn/^ 



Saed 



Samadi 



1950 Lincoln Ave,Apt.605 Montreal QC 
H3H2N8,CANADA 



Number 2 of 



WARNING: Information on this form may become public. Credit card information sliould not 
be Included on this form. Provide credit card Information and authorization on PTO-2038. 



Copy provided bv USPTO from the PACR Imaae Database on 11/21/2003 



Algorithm for Exact Computation of Coefficients of Universal Maximally Flat FIR 

Filters 

Saed Samadi and Aldnori Nishihara 
Abstraf:t 

The universal masdmally flat FIR filters have an exact closcd-from ejq>resmon for the value of their impulse 
response coefficients. The exptesaon involves multiple summattons and binomial coeffidents. The binomial co- 
effidents may posably have non-integ»r arguments. The direct computation of the impulse response coeffidents 
of the ffiters for a given set of design parameters is computationally expeasive and requires espedal software 
routined for evaluation of large binomial coeffidents. We show that the impulse response coefficients of the uni- 
versal maximally flat FIR filters may be compute effidtotfe^ uSinff a twi>-stk|^ algofithin. Fbr a filter of order 
N with K zeros at a = -1 , the fiist stage of the aigorithui uses a two^term recurrence to generate a sequence 
of length N - K + 1 . This sequence undergoes an N-step iterative process that involves additions, subtractions 
«iil diyifflons hy 2. Th^ algonVhin dees n6t involve, any form of binoiial cdbffidents. Moreover, the algorithm 
requires less number of arithmetical operations compared to a dired: evaluation using the closed-form formula. 
The algorithm is mostly suitable for real-tfane deration of the coeffidents on a DSP drip because it may be 
easily implemented to a computational envinwiment that has restrictions on the available hardware or software 
resources. 
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dfAay systems, halfband filters, Bemstdn polyn<»nials, forward shift operator, binomial coefiteJents. 
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I. Introduction 

More thao thirty years have passed since the pubUcation of the seminal paper of O. Hcmnann 
on the design of maximaUy flat FIR digitai filters [1]. Of the body of the related research that has 
been published since Herrmann's article [2)-[211 S the most significant contribution to the theory of 
lowpass masdmaUy flat FIR filters has beraa Bahor's results on tiie design of filters with amultancous 
conditions on ampUtudc and group-dday response [8]. In recent years, it has been recognized that 
Bahcr's filters form a unifying dass of maadmally flat filters. This class encompasses various types 
of seemingly distinct systems [20]. Calling this dass of filters the universal maximaUy flat digital 
filters, the authors simplified Bahor's formula for the transfiar function and showed that lowpass filters 
of even or odd lengths, with Unoar or nonlinear phase response characteristics, as Wdl as fractional 
delay and half-band filters may be readily obtained by setting the do»gn parameters in an appropriate 
manner. Other results show that Baher's universal maadmally fiat filtrars are exact and optimal with 
respect to polynomial signals [21]. Apart fixjm their importance because of this remarkable properly, 
some authors have envisaged other ^pUcation scenarios for maximally flat filters. For instance, [10] 
discusses a situation where the use of an Loo- or Li-optimal filter may not result in a desirable design. 

Notwithstandmg universal maximaUy flat digital filters have an explidt formula for their transfer 
'functions and, thus, may be designed readUy using a computer algebra system, the simplified formulas 
provided in [20] arc still unwiddy. A dkect use of the formulas involves computation of binomial 
coeffidents with possibfy non-inte©sr kidioes that appear as terms of the summand in a three-fold 
sum. Binomial coeffidents with non-intogcr indices cannot be computed by the simple method of 
Pascal's triangle. Another common problem with computations involving the binomial coeffidents is 
due to thdr tendency to introduce very largo integers to the intermediate steps of the computation. 

' We do not claim that Ihe list of references given here b a comprdiensive one. The ciusd papers are those with either 
tiieoteUcal or historical implications. 
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With the computer power al our disposal today, this may not be a serious impediment in many 
situations. However, we can stiU think of cases where Umitations on hardware and softwaxo resources 
dictate the use of efficient means of computation of filter coefficients. It is desirable that such methods 
be free of binomial coefficients and have a simple iterative structure. Moreover, we know that if 
the design parameter associated with the group delay is limited to rational numbers, as is the case 
with many realistic settings, the values of the impulse response coefficients are rational numbers. In 
such cases, it is important to ensure that the answer to the numerical problem of computation of an 
impulse response coefficient is exactly a fraction, say 1/3, not a floating number that may get printed 
as "0.333333335". Therefore, it is desirable for a computation algorithm to maintain this property of 
the impulse response coefficients and yield rational values for rational design parameters.- 

The contribution of this paper in connection with the above computational issues is as follows. To 
avoid a direct use of the computationally ejtpcnsivc formula developed in [20] for the task of cxaict 
computation of the coefficients of universal maximaUy flat filters, we devise a recursive algorithm that 
is free from binomial coefficients of any form, requires less multiplications, and may be implemented 
using rational arithmetic. We start by deriving a simple recurrence relation for generation of a sequence 
that forms the coefficients of a Bemstcin-form representation of the transfer function. The Bcmstem- 
-form representation of the transfer function was fully developed in [20]. Then, a simple method 
for converting the Bcmstem-form polynomials to the power form is introduced. A variance of this 
method is derived to convert the representation of the transfer function given in [20] to the power- 
form representation using successive additions, subtractions and divisions by 2. An algorithm for 
computation of the impulse response coefficients is obtMned by combining our recurrence relation and 
Bemstein-to-power conversion method. 

The treatment is self-contained featuring complete proofs and supported by a diagram and an 
example. The mathematical results required for derivation of the algorithm arc developed m Section 11. 



rtnrwt nrnviffprf hv IIRPTn ffrr%m thn PACR Imsuta risatnhsMA nn i1/9i/9nn!) 



fi.o«-!i--Hi.iB"^:i 3 ,.10.1 ?iiOEi: 



The detailed description of the algorithm together with a pscudo-codc for its implementation arc 
given in Section HI. A concrete example is also worked out in this section. Conclusions arc drawn in 
Section IV where the computational complexity of the algorithm is studied as well. 

II. Mathematical Underpinning 
It was shown in [201 that the transfer function of the universal maximaUy flat digital filters, a family 
of filters identical to those proposed by Bahcr \S], is ^von Ijy 

0<j<N-lC 

where the coefficients bj arc characterized by the threo-tcrpi recurrence 

j.bj + 2dbH,-(j-'N-2)b}-2 = 0. i^l. (2) 

with the initial values 

bo = 1 and b-i = 0. (3) 
By converting (1) to the power form representation 

0<k<N 

•the impulse response coefficients Hk become [20] 

_ ^ ^ (-ir-^(V)(t.'')a)0 N (5) 

0£)<N-K 0<p<1cOSl<J 

Obviously. (5) is a computationally expensive formula. To ameUoratc this situation, wo devise an 
algorithm that exploits tho combination of (1) and (2) to compute hk for Ic = 0.. . . ,N. 

First, let us write (1) in the traditional Bernstein form by mtroducuig the sequence b,' defined by 

I 0<j<N. 
0, otherwise. 
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Then wo have 

This is the Bcmstoin fonn of the transfer function in the traditional sense Although (7) has an 
additional binomial coefficient compaicd to (1), wo will soo later that tho introduction of b/ results in 
considerable simplification. 

A recurrence for b/ is obtained by substituting (6) into (2) and dividing both sides by K?) ?^ 0- 

Thus, 

b; + 2d^b,'_,-(j-N-2)^b(_, = 0. l<j<N. (8) 
The above recurrence can be simplified hy invoking the definition of the binomial coefficients [22] and 
a few simple algebraic steps. It £9}low$ that 

with the initial values 

bS = l and bi, =0. (10) 

Thus, for given values of N and d wc can compute the values of b,' recursively by (9). 

The next step is to convert (7) to the power form (4) and obtain the impulse response coefficients. 
A direct cxpanaon of the summand In (7) usfaig the binomial theorem will introduce additional sums 
abd new binomial coefficients. Hence, this is not a desirable action to take. The other alternative is to 
use a conversion identity weU-known in the mathematics literature and in the field of computer-aided 
geometric design. The identity is based on takmg successive differences of the Bernstein coefficients. 
We state this identity as a proposition and provide a proof for it. 

Pmposttwn i: Let the Bemstem form of the polynormal p(x) = Eo<i<N Pv*^ S^<^ 

P(x)= L bt(^)xMl-xr-v. (11) 

0£i<N \ / 
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Then 

Pi=^^^A^bo. i = 0,....N. (12) 

where the operator A is the forward difference operator. 

The forward diifcrcnco operator is defined hy the gisncral recurrence 

togotiicr with 

A'b» = Abx = bi+i-bi, A«bi = K. (14) 

For instance, Abo = b, - bo and A^bo = bz - 2b, + bo- In this notation, the operator acts on the 
indices, written as subscripts, of the members of a soqucncij. To prove the propoaitionj wc also need 
the forward shift operator E. The effect of E on the indcJt I bf a Sequence like bi is specified by 

Eibi=:E(Ei-ibO. (15) 



together with 



E is related to A by 



E»bi = Eb^ = bt+,, EObt = bi. (16) 



A = E-1. (17) 
Pnof of proposition. We start with the Bernstein form of p(x) and expand the term (1 -x)*^"* using 



the binomial theorem to obtain 

(18) 



By the symmetry of the binomial coefficients and the trinomial revision property {22], wo have 



On the other hand 

bi = E^bo. (20) 

Thus, we can write 

,(x.= E Z i-')'<-'(^)LlTi,)'''-'^''- 

0<i<N 0<j<N-l \ ' / \ / 

The tangos for the indices of the two-fold sum above may be interchanged as 



0<1<N 0<5<N-l 0<j<N 0<i<i 



wtth no effect on the value of the sum. Now, wc substitute N - j for j in the sum and obtain 

0<N-i<NO<t<) \ '/ A ' . _ 



OSN-}<N 0<t<J 

The ri^t side above simpUfies to 



Note that 



It then follows that 



o<i<} ^' ^' 



(22) 



(23) 



(25) 



(26) 



0<)<N V ' / 

that is the power form representation of p(x). Rom the uniqueness of the power from representation 
of a polynomial and the relation between the operators E and A it follows that Vi = (i)A^bo. ■ 

We now present our main and last result in this preparatory phase before proceeding to the derivation 
of the algorithm in the next section. 

Corollary 1: The transfer function HN,K.d(z) « ffiwen by 

HN.K.a(.) = (4^ + ^--^)^t.S. (27) 



Proof. Let 
Then, wo can write 

>N-J. (29) 



H^..K..(x)= H mC;')-'(1-x)^ 



0<j<N-K 

By the above proposition wc have 



0<j<N \^ / 

By substituting for x in tcnns of x"^ in tho above relation and expanding the sununand using the 
binonual theorem, we obtain 

a<i<N \ ' / i-<j \ / 



o<i<N 

By the trinomiai rcviaon identity wc have 

(^)(0=C)(-^)- 

We can also modify the sum indices and write 

0<1<N \ / ^ 

. Using p = j - i as an index for the second sum, we obtain 



Hn 



This rampiifics to 



0<1<N 

By mvoking tho buiomial theorem, we can further write 



2 2 

By ushig the relation between A and E, (27) foUowa 



(31) 



(32) 



(33) 



(34) 



0<1<N \ / 



, > ft ■ ^ ^,-i^Nv (36) 
HN.K.d(z) = (l+ 2 ~ 2* ' 



I 



III. The Algorithm 

The algorithm consists of two major stages. In the first stage, we run the recurrence (9) from j = 1 
uptO| = N-'Kto obtain the sequence 

B' = (1,b{...,b{;,^K»0....,0). (37) 

The entries bj of B' for N— K < j < N arc set to zero. The next step is to expand the transfer function 
using (27). This can be done by noting that 

(!±l + lzl.-.)%J = (l+l + l^.-)((l±l + l^.-)'-«) (38) 

for arbitrary integer p > 0, Eyvaluation of both sides above in the power form yields 

where h|^^ is a sequence of coefficients defined by 

h[P» = (l4^1)^(h«^).h5>....), (40) 
and operator £ effects a forward shift on the index j. Thus, we obtain 

2:hi^»z-^ = £ (^^^"'^ + ^i"'^ ^"^'')- 

It follows that 

h«p) = 1+1 hr'» + 1^ ht-i'\ (42) 

The ranges for i and p arc specified by 

1 < P < 0 < i < p. (43) 

The initial values are given by 

hg>^=B'. h?^ = (0,0....). (44) 
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This completely specifics the rccunencc wc use in the second stage of our algorithm. After completing 
the Nth step, at the point when p haa reached N . the value of the vth impulse response coefficient is 
g^ven by 

hi = first term of sequence h[^\ i=0,...,N. (45) 

The recurrence is illustrated in Fig. 1 for N = 3 and an unspecified value of K, It can be visually verified 
that the recurrence has a triangular structure. A pseudo-code of a procedure for implementation of 
the algorithm is given in Fig. 2. 



definitions 

h 1-D array for impulse response coefficients 
B' 1-D array for bj 

r 3-D array for the intermediate values 
i counter for t as depicted in Fig. 1 
i counter for t and B' 
p raunter for T as depicted in Fig. 1 
procedure 

for 3 = 0 to N step +1 do 
^0 // Initialization 

endfor 

for 1, j,p =s 0 to N step 1 do 

r(p,t,}) ^ 0 // Initialization 
endfor 

B'[ll*-1 // Imtial values 
for toN-Xstep+l do // (A) 
• B'Q] 4- ((2d) B'0 - 11 + (j - 1 ) B'B - 21) 
endfor // END OF (A) 
for j = 0 to N ~ K step 1 do 

t[0,0J1 4- B'Bl // Initial values 
endfor 

for p = 1 to N step 1 do // (B) 
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for i from 0 to p step 1 do 

for j from 0 to N -p step 1 do 

Tlp.iJ]:={r[p-l.i-lJ]-T[p-l.i-lJ + l])/2+(T[p-1,i>i] + T^^ 

endfor 
endfor 

endfor // END OF (B) 
for i = 0 to N step 1 do 
Hli] (-T[N,i.Ol end do 
return H 
endprocedure 

Fig. 2. Pseud<«ode of a procedure that returns tKe values of impulse response coefficients as an 4rray. . 
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The foUowing example helps elucidate the details. Let N = 3, K = 1 and d = (-l)/4. The actual 
values of the sequences i = 0, .... 3, corresponding to the 3-D array r in the procedure of Fig. 2, 
are as follows 

h?^ hi^J 14^^ 14^^^ 

p = 0: (1,1/6,-11/24,0) 

p = l: (7/12,-7/48, ~11/4S) (5/12.5/16,-11/48) (-^6) 
p=2: (7/32.-3/16) (35/48,1/12) (5/96,13/48) 

p = 3: (1/64) (39/64) (31/64) (-7/64) 

Only non-zero values of the sequences i = 1, ... ,3, arc presented. Note that all computations 
have been done in the field of rational numbers in an exact manner. Hence, round-off errors are 
avoided. 

IV. Conclusion and Remarks 

We have shown that an alternative representation for the transfer function of universal maximally 
flat PIR filters qnabies us to dcvcaop a very simple, ajigorithm for evaluation of the impulse response 
coefficients of the filters. The alternative form of the transfer function is based on the shift operator 
E and is derived by means of a technique referred to as finite calculus in [22] . 

The algorithm consists of two stages. In the first stage wo run a two-term recurrence relation and 
dbtain a 1-D sequence of numbers. If the value of the group delay parameter d is a rational munber, 
then the outputs of the recurrence are rational. The outputs are then used in a simple N-step algorithm 
of triangular structure. In each step of this algorithm, the values of tho sequences obtained in the 
preceding step undergo simple additions, subtractions and divisions by 2 to yield sequences that are 
shorter in size. Upon tho completion of the Nth step, N sequences arc obtained. The first entries of 
these sequences are tho values of the impulse response coefficients. As long as the computations in the 
second stage are done in the field of rational numbers the final values remain rational. This provides 
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an accurate and fast means for exact computation of the impulse response coefficients. 

A v«>r8^case analysis of the computational complexity of the proposed algorithm is the final topic 
of this section. The computational complexity be analyzed hy evaluatmg the comple:dty of the 
simple forwloop marked by (A) and the nested for-loop designated by (B) in Fig. 2. For the number 

of additions and subtractions C±, wc obtain 

c^= E ^+ H H L ^- ^^'^ 

l<j<N-K l<p<NO<i<pO<j<N-p 



For a closed-form expression, note that 



that can be evaluated as 



1<p<N0<l<p+l 



C±=4(N-K)+3 ^ (p+2)(N-p). 

0<p<N 



For the number of multiplications C*, we have 

C*= ^ 2 = 2(N-K), 

|<J<N-K 

' The number of divisions C/ is given by 



(48) 



(49) 



Thus, wo get 

C±=4(N-X) + 3(|N + N2 + iN3). (50) 



(51) 



C/= L 1+ H L H ^ ^'"^ 

1<JSN-K lSP<NOSl<pO<}<N-n> 



(53) 



that becomes 

C/ = (N-K)+2(|n + N* + ^N3). 
In comparison, for the direct evaluation using (5). the number of additions and subtractions cg"^* 

is 

cDiroct> ^ L Z ^ 
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The right side above is of the order O(N-). The number of rnultipUcaiions C?^', on the other 
hand, is 

Cpirect> Y_ H L'^- 

0<i<N-K 0<p<kO<l<J 

The right ^de above is of the order OIN^). Fbr the nun*er of divisions in a direct evaluation. ^ can 
write 

Cpirect> £ H £4. (56) 

0<i<N-K 0<p<lc 0<i<i 

The right side above, again, is of order O(N^). Thus, in addition to the simplicily of most of the 
division operations in thi. proposed algprithm (that axe divisions by 2). a significant reduction of the 
overaU complffldty is accommidatod in terms of sheer numbers. 
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Fig. 1. Illustrative representation of the algorithm. Each node represents a sequence It^^ . The nodes marked 
by 0*8 represent zero sequences. 



